Bienvenides a la cuarta tarea del curso Statistical Thinking. Esta tarea tiene como objetivo evaluar los contenidos teóricos de la primera parte del curso, los cuales se enfocan principalmente en introducirlos en la estadística bayesiana. Si aún no han visto las clases, se recomienda visitar los enlaces de las referencias.
La tarea consta de una parte práctica con el fin de introducirlos a la programación en R enfocada en el análisis estadístico de datos.
Slides de las clases:
Videos de las clases:
Documentación:
En la siguiente sección deberá resolver cada uno de los experimentos computacionales a través de la programación en R. Para esto se le aconseja que cree funciones en R, ya que le facilitará la ejecución de gran parte de lo solicitado.
Para el desarrollo preste mucha atención en los enunciados, ya que se le solicitará la implementación de métodos sin uso de funciones predefinidas. Por otro lado, Las librerías permitidas para desarrollar de la tarea 4 son las siguientes:
# Manipulación y funcional
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.2 ✔ tibble 3.3.0
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.1.0
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(dplyr)
library(tidyr)
library(purrr)
# Plots
library(scatterplot3d)
## Warning: package 'scatterplot3d' was built under R version 4.5.2
library(ggplot2)
library(plotly)
##
## Attaching package: 'plotly'
##
## The following object is masked from 'package:ggplot2':
##
## last_plot
##
## The following object is masked from 'package:stats':
##
## filter
##
## The following object is masked from 'package:graphics':
##
## layout
# Varios plots
library(gridExtra)
##
## Attaching package: 'gridExtra'
##
## The following object is masked from 'package:dplyr':
##
## combine
# Análisis bayesiano
library(rethinking)
## Loading required package: rstan
## Warning: package 'rstan' was built under R version 4.5.2
## Loading required package: StanHeaders
## Warning: package 'StanHeaders' was built under R version 4.5.2
##
## rstan version 2.32.7 (Stan version 2.32.2)
##
## For execution on a local, multicore CPU with excess RAM we recommend calling
## options(mc.cores = parallel::detectCores()).
## To avoid recompilation of unchanged Stan programs, we recommend calling
## rstan_options(auto_write = TRUE)
## For within-chain threading using `reduce_sum()` or `map_rect()` Stan functions,
## change `threads_per_chain` option:
## rstan_options(threads_per_chain = 1)
##
## Do not specify '-march=native' in 'LOCAL_CPPFLAGS' or a Makevars file
##
## Attaching package: 'rstan'
##
## The following object is masked from 'package:tidyr':
##
## extract
##
## Loading required package: parallel
## Loading required package: dagitty
## rethinking (Version 2.01)
##
## Attaching package: 'rethinking'
##
## The following object is masked from 'package:purrr':
##
## map
##
## The following object is masked from 'package:stats':
##
## rstudent
Si no tiene instalada la librería “rethinking”, ejecute las siguientes líneas de código para instalar la librería:
install.packages("rethinking")
En caso de tener problemas al momento de instalar la librería con el código anterior, utilice las siguiente chunk:
install.packages(c("mvtnorm","loo","coda"), repos="https://cloud.r-project.org/",dependencies=TRUE)
options(repos=c(getOption('repos'), rethinking='http://xcelab.net/R'))
install.packages('rethinking',type='source')
Por último, si sigue con problemas, intentar con el siguiente script:
https://github.com/dccuchile/CC6104/blob/master/scripts/rethinking_install.R
En caso de no lograr instalar la librería a pesar de estos pasos, favor contactar al auxiliar por correo.
Se dispone del archivo local
online_shoppers_intention-2.csv, donde la variable
Revenue indica si un usuario realizó una compra
(TRUE/FALSE o 1/0). Se busca estimar la
tasa de conversión promedio y analizar cómo cambia la
inferencia bayesiana al variar el prior y el tamaño de muestra.
Se pide:
Revenue a variable
binaria {0,1}. Calcule la tasa global de conversión.
quap) bajo prior uniforme y compare sus resultados con el
método de grilla.
¿Son similares estos intervalos? ¿Porqué?
¿Que diferencia tiene con respecto al intervalo de confianza de el enfoque frecuentista?
Respuesta:
Revenue a variable
binaria {0,1}. Calcule la tasa global de conversión.
#Fijamos la semilla para mantener los mismos resultados
set.seed(3)
datos_ecom <- read_csv("datos/online_shoppers_intention-2.csv", show_col_types = FALSE)
datos_ecom$Revenue <- ifelse(datos_ecom$Revenue == "FALSE", 0, 1)
suma <- sum(datos_ecom$Revenue)
total <- nrow(datos_ecom)
prop <- suma/total
#Así está bien o en porcentaje?
cat("Con un total de ", total, " datos, y una cantidad ", suma, " de compras, ")
## Con un total de 12330 datos, y una cantidad 1908 de compras,
cat("la proporción de usuarios que sí realizó la compra es de: ", round(prop*100, 2), "%")
## la proporción de usuarios que sí realizó la compra es de: 15.47 %
Prior A: \(p \sim N(0.14, 0.05)\)
Prior B: \(p \sim \text{Uniforme}(0,1)\).
#Definimos la grid, que son los posibles valores de p
p_grid <- seq(from=0, to=1, length.out=1000)
#Definimos los priors
prior_A <- dnorm(p_grid, 0.14, 0.05)
prior_B <- rep(1, 1000)¿Cómo cambia la forma de la posterior al aumentar \(n\)?
#Los tamaños de muestras que nos piden
muestras <- c(10, 50, 500, 1000)
posterior <- function(prior, sec, nombre){
par(mfrow=c(2,2))
for (n in sec){
muestra <- sample(datos_ecom$Revenue, size = n, replace = FALSE)
x <- sum(muestra)
#probabilidad de observar exactamente x éxitos en n muestras para cada una de las p probabilidades de p_grid
likelihood <- dbinom(x, n, prob=p_grid)
#posterior sin estandarizar
unstd.posterior <- likelihood * prior
#la estandarizamos y tenemos la distribución actualizada de p dada nuestros datos
posteriors <- unstd.posterior / sum(unstd.posterior)
max_post <- max(posteriors)
max_index <- which.max(posteriors)
max_p <- p_grid[max_index]
print(paste("Máximo de posterior para ", nombre, " con n = ", n, ": ", round(max_post, 4), " en p = ", round(max_p, 4)))
plot(p_grid, posteriors,
xlab="Probabilidad de realizar la compra",
ylab="Probabilidad posterior",
main = paste("Posterior con prior", nombre , "- n =", n),
xlim = c(0, 1), ylim = c(0, max(posteriors)))
abline(v=prop, col="red", lwd=2, lty=2)
legend("topright", legend=c("Posterior", "Proporción"),
col=c("black","red"), lty=c(1,2), lwd=2)
}
}
posterior(prior_A, muestras, "A")
## [1] "Máximo de posterior para A con n = 10 : 0.0088 en p = 0.1331"
## [1] "Máximo de posterior para A con n = 50 : 0.0115 en p = 0.1502"
## [1] "Máximo de posterior para A con n = 500 : 0.0265 en p = 0.1471"
## [1] "Máximo de posterior para A con n = 1000 : 0.0346 en p = 0.1712"
posterior(prior_B, muestras, "B")
## [1] "Máximo de posterior para B con n = 10 : 0.0029 en p = 0.3003"
## [1] "Máximo de posterior para B con n = 50 : 0.0104 en p = 0.0801"
## [1] "Máximo de posterior para B con n = 500 : 0.0246 en p = 0.1562"
## [1] "Máximo de posterior para B con n = 1000 : 0.0348 en p = 0.1562"
¿Cómo cambia la forma de la posterior al aumentar \(n\)?
Al ir aumentando n, la posterior se va pareciendo cada vez más a la posterior real. Esto se puede observar al ver que el máximo del gráfico se va acercando cada vez más a 0.15, la proporción que sacamos anteriormente. Esto se puede explicar porque a medida que aumenta la cantidad de datos muestreados, los datos predominan por sobre los priors, concentrando la curva alrededor de la proporción real.
Compararando ambos priors, el prior A varía menos a medida que aumenta n, convergiendo suavemente a la distribución real. Esto se puede explicar, dado que parte con una media en 0.14 muy cercana a la proporción real de unos 0.15. Mientras que el prior B, al inicio tiene menos certeza (todos los valores son igual de probables), entonces con menos datos tiene menos información para acercarse al valor real, pero a medida que aumenta n, va pareciéndose más a dicha proporción.
quap) bajo prior uniforme y compare sus resultados con el
método de grilla.#Para poder comparar los cuatro gráficos de una
par(mfrow=c(2,2))
for (n in muestras){
muestra <- sample(datos_ecom$Revenue, size = n, replace = TRUE)
x <- sum(muestra)
#Como lo vimos en clase
globe.qa <- quap(
alist(
x ~ dbinom(n, p), # binomial likelihood
#Por problemas con el knit cambiamos el dunif(0,1) por dbeta(2,2), que es lo más cercano a uniforme posible que no tenga problemas con x = 0 o x = n
p ~ dbeta(2,2) # uniform prior
) ,
data=list(x=x, n=n))
# Extraemos el valor de p desde el globe.qa
posterior_quap <- extract.samples(globe.qa)
#Densidad de p obtenida de la aproximación
densidad<- density(posterior_quap$p)
max_densidad <- max(densidad$y)
max_p <- densidad$x[which.max(densidad$y)]
print(paste("Máxima densidad de Laplace: ", round(max_densidad, 4), " en p = ", round(max_p, 4)))
plot(densidad, main = paste("Posterior con Laplace - n =", n),
xlab = "Probabilidad de compra", ylab = "Densidad",
xlim = c(0, 1), ylim = c(0, max(densidad$y)))
abline(v=prop, col="red", lwd=2, lty=2)
legend("topright", legend=c("Posterior", "Proporción"),
col=c("black","red"), lty=c(1,2), lwd=2)
}
## [1] "Máxima densidad de Laplace: 5.1134 en p = 0.0762"
## [1] "Máxima densidad de Laplace: 10.812 en p = 0.0832"
## [1] "Máxima densidad de Laplace: 25.0051 en p = 0.1424"
## [1] "Máxima densidad de Laplace: 34.8464 en p = 0.1527"
Al comparar Laplace con la grid approximation, ambos métodos producen resultados coherentes en general. Se observan algunas diferencias cuando n es pequeño y podría deberse principalmente a la variabilidad de la muestra, porque estamos analizando muestras aleatorias distintas. A pesar de esto, con n grandes (500 y 1000) la posterior se concentra y la aproximación de Laplace coincide muy bien con la solución por grid.
Veamos cómo se vería si comparamos con muestras iguales
comparacion <- function(){
for (n in muestras){
cat("Para n: ", n)
# Usamos la misma muestra
muestra <- sample(datos_ecom$Revenue, size = n, replace = FALSE)
x <- sum(muestra)
# grid
likelihood <- dbinom(x, n, prob = p_grid)
unstd_post <- likelihood * prior_B
posteriors <- unstd_post / sum(unstd_post)
# MAP grid
MAP_grid <- p_grid[which.max(posteriors)]
cat("MAP (Grid): ", round(MAP_grid, 4), "\n")
#laplace
modelo <- quap(
alist(
x ~ dbinom(n, p),
p ~ dunif(0,1)
),
data = list(x = x, n = n)
)
muestras_quap <- extract.samples(modelo)
dens_quap <- density(muestras_quap$p)
# MAP
MAP_laplace <- dens_quap$x[which.max(dens_quap$y)]
cat("MAP (Laplace): ", round(MAP_laplace, 4), "\n")
plot(
p_grid, posteriors, type = "l", col = "blue",
xlab = "p", ylab = "Posterior",
main = paste("Grid vs Laplace (n =", n, ")")
)
lines(dens_quap$x, dens_quap$y / max(dens_quap$y) * max(posteriors),
col = "red", lwd = 2)
legend("topright",
legend = c("Grid", "Laplace"),
col = c("blue", "red"), lwd = 2)
}
}
comparacion()
## Para n: 10MAP (Grid): 0.3003
## MAP (Laplace): 0.3044
## Para n: 50MAP (Grid): 0.1201
## MAP (Laplace): 0.1206
## Para n: 500MAP (Grid): 0.1481
## MAP (Laplace): 0.1514
## Para n: 1000MAP (Grid): 0.1572
## MAP (Laplace): 0.158
Podemos observar que ahora ambas técnicas producen resultados aún más
coherentes, donde los MAP son parecidos entre ambas y se parecen aún más
cuando n aumenta.
¿Son similares estos intervalos? ¿Porqué?
¿Que diferencia tiene con respecto al intervalo de confianza de el enfoque frecuentista?
par(mfrow=c(2,2))
for (n in muestras){
muestra <- sample(datos_ecom$Revenue, size = n, replace = TRUE)
x <- sum(muestra)
globe.qa <- quap(
alist(
x ~ dbinom(n, p),
p ~ dunif(0, 1)
),
data=list(x=x, n=n)
)
posterior_quap <- extract.samples(globe.qa)
cred_interval_50 <- PI(posterior_quap$p, prob=0.5)
cred_interval_95 <- PI(posterior_quap$p, prob=0.95)
hpdi_50 <- HPDI(posterior_quap$p, prob=0.5)
hpdi_95 <- HPDI(posterior_quap$p, prob=0.95)
cat("Intervalo de Credibilidad 50%: ", cred_interval_50, "\n")
cat("Intervalo de Credibilidad 95%: ", cred_interval_95, "\n")
cat("HPDI 50%: ", hpdi_50, "\n")
cat("HPDI 95%: ", hpdi_95, "\n")
densidad <- density(posterior_quap$p)
plot(densidad, main = paste("Posterior con Laplace - n =", n),
xlab = "Probabilidad de compra", ylab = "Densidad",
xlim = c(0, 1), ylim = c(0, max(densidad$y)))
abline(v = cred_interval_50, col = "red", lwd = 2, lty = 2)
abline(v = cred_interval_95, col = "blue", lwd = 2, lty = 2)
abline(v = hpdi_50, col = "green", lwd = 2) # HPDI 50%
abline(v = hpdi_95, col = "purple", lwd = 2) # HPDI 95%
legend("topright",
legend = c("PI 50%", "PI 95%", "HPDI 50%", "HPDI 95%"),
col = c("red", "blue", "green", "purple"),
lty = c(2, 2, 1, 1),
lwd = 2,
bty = "n",
cex=0.7)
}
## Intervalo de Credibilidad 50%: 0.1166846 0.284856
## Intervalo de Credibilidad 95%: -0.04031261 0.4502999
## HPDI 50%: 0.1159058 0.2837031
## HPDI 95%: -0.04045292 0.4493685
## Intervalo de Credibilidad 50%: 0.08907396 0.1514695
## Intervalo de Credibilidad 95%: 0.0306897 0.2108611
## HPDI 50%: 0.09439276 0.1562724
## HPDI 95%: 0.02904603 0.2086294
## Intervalo de Credibilidad 50%: 0.1376426 0.1588014
## Intervalo de Credibilidad 95%: 0.1165092 0.1788537
## HPDI 50%: 0.1376402 0.1587902
## HPDI 95%: 0.1177929 0.1798519
## Intervalo de Credibilidad 50%: 0.1444415 0.1596034
## Intervalo de Credibilidad 95%: 0.1298629 0.1746387
## HPDI 50%: 0.1438393 0.1589521
## HPDI 95%: 0.1299458 0.1747096
¿Son similares estos intervalos? ¿Porqué?
Sí, los intervalos son similares. Se puede observar que en algunos casos, los colores de los PI 50% y 95% casi ni se ven porque son tapados por los de HPDI 50% y 95% correspondientes. Esto ocurre porque la forma de la posterior es bastante simétrica y unimodal. Entonces al calcular el intervalo con ambos métodos no se encuentran mayores diferencias, puesto que el intervalo más estrecho está alrededor del centro de forma más o menos simétrica y así el método que deja colas iguales y el que busca el intervalo más estrecho son bastante similares, y están ambos centrados alrededor de la moda.
¿Que diferencia tiene con respecto al intervalo de confianza de el enfoque frecuentista?
El intervalo de confianza en el enfoque frecuentista indica los límites en los cuáles se encontrará el valor buscado un x% de las veces, si es que el experimento se realiza infinitas veces. Mientras que los intervalos de credibilidad reportan un intervalo que contiene la creencia de encontrar con un x% de probabilidad el valor real de ese parámetro ahí.
A work by CC6104